function MCMC!(theta,dataset,counterfactuals,priors,iter::AbstractRange;
			save_params = minimum(iter):10:maximum(iter),
			save_lv = minimum(iter):100:maximum(iter),
			save_lv_all = [maximum(iter)],
			save_counterfactuals=save_lv,
			outputpath=outputpath, history = typeof(theta)[],
			callback = ((theta,dataset)->false), model=model)
	for year in keys(dataset)
		(cc,lv,temps) = dataset[year]
		precomputeMeanUtilities!(theta,lv,cc,temps)
	end
	for t in iter
		println("iter $t")
		iterate!(theta,dataset,priors,model)
		if t in save_params
			println("saving params")
			push!(history,deepcopy(theta))
			params_out = Dict{String,Any}(string(k)=>v for (k,v) in pairs(theta))
			matwrite(outputpath*"/params_$t.mat",params_out)
		end
		if t in save_lv
			println("saving latent variables")
			if t in save_lv_all
				yrs = [2010,2011,2012]
			else
				yrs = [2012,]
			end
			for year in yrs
				(_,lv,_) = dataset[year]
				lv_out = lv
				JLD2.@save(outputpath*"/latentvariables$(year)_$t.mat",lv_out)
			end
		end
		if t in save_counterfactuals
			println("running counterfactuals")
			ww = runCounterfactuals(theta,dataset,counterfactuals)
			CSV.write(outputpath*"/counterfactuals_$t.csv",ww)
		end
		callback(theta,dataset) && break
	end
	history
end

########################################################
## single iteration
########################################################
function iterate!(theta,dataset,priors,model)
	ntypes = dataset[2011][1].ntypes
	for year in keys(dataset)
		(cc,lv,temps) = dataset[year]
		(J,I) = size(lv.U)
		Threads.@threads for ii = 1:I
			updatemu!(theta,lv,cc,temps,ii)
	        updateIndividual!(theta,lv,cc,temps,ii) #(1) utility and availability
			updateRC!(theta,lv,cc,temps,ii) #(2) students' random coefficients on program char.
	    end
	end
	for typ=1:ntypes
		updateVCVofRCs!(theta,dataset,typ,priors.VCVrc)
		updateOutsideOptionParameters!(theta,dataset,typ,priors)
		updateProgramRE!(theta,dataset,typ,model) #(a)  program effects
		updateFixedTerms!(theta,dataset,typ,priors.fixed)	#(b) coef on j (fixed) terms
		updateMatchTerms!(theta,dataset,typ,priors.ij)	#(c) coef on ij terms
		updateVarianceOfProgramRE!(theta,dataset,typ,model)  #(d) variance of program random effects (if RE model)
		updateAlpha!(theta,dataset,typ,priors.alpha) #(5) pr of unavailability
		updateOutcomeTerms!(theta,dataset,typ,priors.outcome)
	end
end

function updateOutsideOptionParameters!(theta,dataset,typ,priors)
	#oo variances
	prior_varU0 = priors.var_oo0 #InverseGamma(1,1)
	prior_varU1 = priors.var_oo1 #InverseGamma(1,1)
	for year in keys(dataset)
		(cc,lv,temps) = dataset[year]
		inds_typ = findall(cc.typ .== typ)
		oo0 = lv.U0[1,inds_typ]
		oo1 = lv.U0[2,inds_typ]
		prior_varU0 = BayesVariance(oo0,temps.mu0[1,inds_typ],prior_varU0)
		prior_varU1 = BayesVariance(oo1,temps.mu0[2,inds_typ],prior_varU1)
	end
	theta.sigsqOO0[typ] = rand(myRNG[Threads.threadid()],prior_varU0)
	theta.sigsqOO1[typ] = rand(myRNG[Threads.threadid()],prior_varU1)
	#oo mean
	prior_muU0 = priors.oo0
	prior_muU1 = priors.oo1
	for year in keys(dataset)
		(cc,lv,temps) = dataset[year]
		inds_typ = findall(cc.typ.==typ)
		oo0 = lv.U0[1,inds_typ]
		oo1 = lv.U0[2,inds_typ]
		Xoo = cc.Xoo[inds_typ,:]
		prior_muU0 = BayesGLS(oo0, theta.sigsqOO1[typ], Xoo, prior_muU0)
		prior_muU1 = BayesGLS(oo1, theta.sigsqOO0[typ], Xoo, prior_muU1)
	end
	theta.betaOO0[:,typ] .= rand(myRNG[Threads.threadid()],prior_muU0)
	theta.betaOO1[:,typ] .= rand(myRNG[Threads.threadid()],prior_muU1)
	#fill in
	for year in keys(dataset)
		(cc,lv,temps) = dataset[year]
		inds_typ = findall(cc.typ.==typ)
		Threads.@threads for ii in inds_typ
			temps.mu0[1,ii] = dot(view(theta.betaOO0,:,typ),view(cc.Xoo,ii,:))
			temps.mu0[2,ii] = dot(view(theta.betaOO1,:,typ),view(cc.Xoo,ii,:))
		end
	end
end

function updateVarianceOfProgramRE!(theta,dataset,typ,::ProgramFE)
	theta.sigsqProgramRE[typ] = 0
end

function updateVarianceOfProgramRE!(theta,dataset,typ,::ProgramRE)
	theta.sigsqProgramRE[typ] = rand(myRNG[Threads.threadid()],BayesVariance(theta.programRE[:,typ],0))
	return
end

function updateAlpha!(theta,dataset,typ,priorAlpha)
	posteriorAlpha = priorAlpha
	for year in keys(dataset)
		(cc,lv,temps) = dataset[year]
		posteriorAlpha = _updateAlphaYear!(cc,lv,temps,posteriorAlpha,typ)
	end
	theta.alpha[:,typ] .= rand(myRNG[Threads.threadid()],posteriorAlpha)
end

function _updateAlphaYear!(cc,lv,temps,posteriorAlpha,typ)
	inds = cc.inds_WL[typ]
	Astar_typ = temps.Astar_typ[typ]
	for (n,ind) in enumerate(inds)
		Astar_typ[n] = lv.Astar[ind]
	end
	 BayesOLS(Astar_typ,cc.Xfriction[typ], posteriorAlpha)
end


#Deterministic function to compute utility components.
function precomputeMeanUtilities!(theta,lv,cc,temps)
	for typ = 1:cc.ntypes
		mul!(view(temps.ufe,:,typ),cc.Xfixed,view(theta.beta_fixed,:,typ)) #cc.Xfixed*theta.delta_Xfixed[:,typ]
		temps.ure[:,typ] .= theta.programRE[cc.programIDs,typ]
	end
	Threads.@threads for ii=1:size(lv.U,2)
		mytyp = cc.typ[ii]
		temps.mu0[1,ii] = dot(view(theta.betaOO0,:,mytyp),view(cc.Xoo,ii,:))
		temps.mu0[2,ii] = dot(view(theta.betaOO1,:,mytyp),view(cc.Xoo,ii,:))
		mul!(view(temps.urc,:,ii),cc.Xrc,view(lv.rc,:,ii))
		mul!(view(temps.umatch,:,ii),cc.Xij[ii],view(theta.beta_ij,:,cc.typ[ii]))
	end
end

#updateRC!: use mui as temp workspace to store u-deterministicpart
#regress (u-deterministicpart) on terms with RC's.
function updateRC!(theta,lv,cc,temps,ii)
	updatemu!(theta,lv,cc,temps,ii,true,true,true,false)
	typ = cc.typ[ii]
	nrc = size(cc.Xrc,2)
	mui = view(temps.mu,:,ii)
	mui .= view(lv.U,:,ii) .- mui
	posterior = BayesOLS(mui,cc.Xrc,MvNormal(zeros(nrc),theta.Sigma_rc[:,:,typ]))
	lv.rc[:,ii] .= rand(myRNG[Threads.threadid()],posterior)
	mul!(view(temps.urc,:,ii),cc.Xrc,view(lv.rc,:,ii))
end

function updateVCVofRCs!(theta,dataset,typ,prior)
	for year in keys(dataset)
		(cc,lv,_) = dataset[year]
		inds = findall(cc.typ .== typ)
		prior = BayesVCV(lv.rc,prior,inds)
	end
	theta.Sigma_rc[:,:,typ] .= rand(myRNG[Threads.threadid()],prior)
end

function updateProgramRE!(theta,dataset,typ,::ProgramFE)
	theta.programRE[:,typ] .= 0
	for year in keys(dataset)
		(cc,lv,temps) = dataset[year]
		temps.ure[:,typ] .= 0
	end
end

function updateProgramRE!(theta,dataset,typ,::ProgramRE)
	nv = 0
	nre = size(theta.programRE,1)
	sumv = zeros(nre)
	nj = zeros(Int,nre)
	for year in keys(dataset)
		(cc,lv,temps) = dataset[year]
		(_nv,_sumv,_nj) = _SufficientStatsProgramRE(lv,cc,temps,typ,nre)
		nv += _nv
		sumv .+= _sumv
		nj .+= _nj
	end
	prior = Normal(0,sqrt(theta.sigsqProgramRE[typ]))
	Threads.@threads for ind=1:nre
		theta.programRE[ind,typ] = rand(myRNG[Threads.threadid()],BayesNormalMean(nv*nj[ind],sumv[ind],1.0,prior))
	end
	for year in keys(dataset)
		(cc,lv,temps) = dataset[year]
		temps.ure[:,typ] .= theta.programRE[cc.programIDs,typ]
	end
end
function _SufficientStatsProgramRE(lv,cc,temps,typ,nre)
	(J,II) = size(lv.U)
	inds = findall(cc.typ .== typ)
	Threads.@threads for ii in inds
		updatemu!(theta,lv,cc,temps,ii,false,true,true,true)
	end
	nv = length(inds)
	sumv = zeros(nre)
	nj = zeros(Int,nre)
	_sumv_j = zeros(J)
	for ii in inds
		@views _sumv_j .+= lv.U[:,ii] .- temps.mu[:,ii]
	end
	for jj=1:J
		ind = cc.programIDs[jj]
		sumv[ind] += _sumv_j[jj]
		nj[ind] += 1
	end
	nv,sumv,nj
end

function updateFixedTerms!(theta,dataset,typ,prior)
	(deltabar,Ainv) = params(prior)
	ndel = length(deltabar)
	XtX = zeros(ndel,ndel)
	XtH = zeros(ndel)
	for year in keys(dataset)
		(cc,lv,temps) = dataset[year]
		_FixedTermsSuffStats!(XtX,XtH,cc,lv,temps,typ)
	end
	A = inv(Ainv)
	V = inv(XtX .+ A)
	deltatil = V*(XtH .+ A*deltabar)
	V_herm = (V + V')/2
	posterior = MvNormal(deltatil,V_herm)
	theta.beta_fixed[:,typ] .= rand(myRNG[Threads.threadid()],posterior)
	for year in keys(dataset) #precalculate "fixed terms" component of utility
		(cc,lv,temps) = dataset[year]
		mul!(view(temps.ufe,:,typ),cc.Xfixed,view(theta.beta_fixed,:,typ))
	end
end

function _FixedTermsSuffStats!(XtX,XtH,cc,lv,temps,typ)
	inds::Array{Int,1} = findall(cc.typ .== typ)
	N = length(inds)
	J = size(lv.U,1)
	Threads.@threads for ii in inds
		updatemu!(theta,lv,cc,temps,ii,true,false,true,true)
	end
	mul!(XtX,cc.Xfixed',cc.Xfixed,N,1)
	mul!(XtH,cc.Xfixed', vec(sum(view(lv.U,:,inds),dims=2) .- sum(view(temps.mu,:,inds),dims=2)),1,1)
end

#these functions are the same as BayesOLS w/ some shortcuts to exploit structure of X's
# function _updateFixedTerms!(theta,lv,cc,temps,typ,prior::MvNormal)
# 	(deltabar,Ainv) = params(prior)
# 	inds = findall(cc.typ .== typ)
# 	N = length(inds)
# 	J = size(lv.U,1)
# 	Threads.@threads for ii in inds
# 		updatemu!(theta,lv,cc,temps,ii,true,false,true,true)
# 	end
# 	A = Array{Float64,2}(inv(Ainv))
# 	V = inv(N .* cc.Xfixed'*cc.Xfixed .+ A) #if A is known, can actually compute this outside the loop
# 	deltatil = V*( cc.Xfixed'*vec(sum(view(lv.U,:,inds),dims=2) .- sum(view(temps.mu,:,inds),dims=2)) .+ A*deltabar) #many terms here can precompute as well
# 	V_herm = (V + V') / 2
# 	MvNormal(deltatil,V_herm)
# end

function updateMatchTerms!(theta,dataset,typ,prior)
	for year in keys(dataset)
		(cc,lv,temps) = dataset[year]
		prior = _updateMatchTerms!(theta,lv,cc,temps,typ,prior)
	end
	theta.beta_ij[:,typ] .= rand(myRNG[Threads.threadid()],prior)
	for year in keys(dataset) #update precalculated "ij terms" utility component
		(cc,lv,temps) = dataset[year]
		Threads.@threads for ii=findall(cc.typ .== typ)
			mul!(view(temps.umatch,:,ii),cc.Xij[ii],view(theta.beta_ij,:,typ))
		end
	end
end
function _updateMatchTerms!(theta,lv,cc,temps,typ,prior::MvNormal)
	(deltabar,Ainv) = params(prior)
	A = Array{Float64,2}(inv(Ainv))
	inds = findall(cc.typ .== typ)
	N = length(inds)
	J = size(lv.U,1)
	ndelta = length(deltabar)
	Threads.@threads for ii in inds
		updatemu!(theta,lv,cc,temps,ii,true,true,false,true)
	end
	Vinv = A + cc.Xij_t_Xij[typ]
	#XU = Vinv*deltabar
	XUs = [zeros(ndelta) for thr=1:Threads.nthreads()]
	uminusmu = [zeros(J) for thr=1:Threads.nthreads()]
	Threads.@threads for ii in inds
		#mul!(Vinv, cc.Xij[ii]', cc.Xij[ii], 1,1)
		#mul!(XU, cc.Xij[ii]', view(lv.U,:,ii),1,1)
		uminusmu[Threads.threadid()] .= view(lv.U,:,ii).-view(temps.mu,:,ii)
		mul!(XUs[Threads.threadid()], cc.Xij[ii]', uminusmu[Threads.threadid()] ,1,1)
	end
	XU = sum(XUs) + A*deltabar
	V = inv(Vinv)
	deltatil = V*XU
	V_herm = (V+V')/2
	MvNormal(deltatil,V_herm)
end

function updateOutcomeTerms!(theta,dataset,typ,prior)
	(deltabar,Ainv) = params(prior)
	ndel = length(deltabar)
	XtX = zeros(ndel,ndel)
	XtH = zeros(ndel)
	for year in keys(dataset)
		(cc,lv,temps) = dataset[year]
		_AddOutcomeSuffStats!(XtX,XtH,cc,lv,temps,typ)
	end
	A = inv(Ainv)
	V = inv(XtX .+ A)
	deltatil = V*(XtH .+ A*deltabar)
	V_herm = (V + V') / 2
	posterior = MvNormal(deltatil,V_herm)
	theta.betaOutcome[:,typ] .= rand(myRNG[Threads.threadid()],posterior)
end

function _AddOutcomeSuffStats!(XtX,XtH,cc,lv,temps,typ)
	inds = findall( (cc.typ .== typ) .& (cc.enroll .> 0) .& .!(cc.specialAdmit))
	X = temps.Xoutcome[typ]
	mul!(XtX,X',X,1,1)
	mul!(XtH,X',view(lv.H,inds),1,1)
end

########################################################
##  regression coefficients
########################################################
"""delta = BayesOLS(A,deltabar,ustar,Xstar):
Draw coefficients from regression of ustar on Xstar.
Assumes that regression errors ustar-Xstar'delta are distributed N(0,1) iid
Prior is delta ~ MVN(deltabar,A^-1).
"""
function BayesOLS(ustar,Xstar,prior::MvNormal)
	(deltabar,Ainv) = params(prior)
	A = inv(Ainv)
	V = inv(Xstar'*Xstar .+ A)  #Xstar has dimensions (I*J,nDelta)
	deltatil = V*(Xstar'*vec(ustar) .+ A*deltabar)  #from AS appendix
	V_herm = (V + V')/2
	posterior = MvNormal(deltatil,V_herm)
end

"""delta = BayesGLS(A,deltabar,u,Sigma,X):
Draw coefficients from regression of u on x where
 u_i = x[ii]'delta + e_i, with e_i ~ N(0,Sigma).
Prior is given by delta ~ MVN(deltabar,A^-1).
"""
function BayesGLS(u,Sigma,X,prior::MvNormal)
	C = getInvCholeskyFactor(Sigma)
	ustar = C'*u
	Xstar = C'*X
	#Xstar = vcat([C'*X[ii] for ii=1:length(X)]...) #(I*J,nDelta)
	BayesOLS(ustar,Xstar,prior)
end

getInvCholeskyFactor(Sigma) = cholesky(Hermitian(inv(Sigma))).U #(J,J) upper triangular cholesky factor such that C'*C==Σ
getInvCholeskyFactor(sigma2::Number) = 1/sqrt(sigma2)

########################################################
##  means and variances
########################################################
"""BayesVariance: posterior variance of normal random variable v with mean mu."""
function BayesVariance(v, mu, prior = InverseGamma(1,1))
	(initial_df,B0) = params(prior)
	posterior = InverseGamma( initial_df+length(v)/2, B0+(sum((v.-mu).^2))/2)
end

"""BayesMean: mean of iid normal random variable v with variance sigma2"""
function BayesNormalMean(v::AbstractArray,sigma2::Number,prior::Normal=Normal(0.0,10.0))
	BayesNormalMean(length(v),sum(v),sigma2,prior)
end
function BayesNormalMean(lengthv::Number,sumv::Number,sigma2::Number,prior::Normal=Normal(0.0,10.0))
	(mu0,sig0) = params(prior)
	term1 = (1/sig0^2 + lengthv/sigma2)
	posterior = Normal( (mu0/sig0^2 + sumv/sigma2)/term1, sqrt(1/term1))
end

"""BayesBinomialMean: mean of iid bernoulli rv's"""
function BayesBinomialMean(a,prior::Beta=Beta(1,1))
	suma = sum(a)
	lengtha = length(a)
	BayesBinomialMean(suma,lengtha,prior)
end
function BayesBinomialMean(suma,lengtha,prior::Beta)
	(α,b) = params(prior)
	posterior = Beta(α + suma, b+lengtha - suma)
end

"""BayesVCV: draw posterior VCV matrix of u~MultivariateNormal with mean 0"""
function BayesVCV(u,prior::InverseWishart=InverseWishart(4,diagm(ones(size(u,1)))),inds=1:size(u,2))
	nRC = size(u,1)
	nInd = length(inds)
	df,Ψ0 = params(prior)
	Ψ = Matrix{Float64}(Ψ0)
	eps_i = zeros(nRC)
	for ii = inds
		@views eps_i .= u[:,ii] #if mean is not zero, would need to subtract mean here!
		mul!(Ψ,eps_i,eps_i',1.0,1.0) #Ψ += eps_i*eps_i'
	end
	posterior = InverseWishart(nInd+df,Ψ)
end
#"""BayesVCV: draw posterior VCV matrix of u~MultivariateNormal with mean 0"""
# function BayesVCV(u,mu,prior::InverseWishart=InverseWishart(4,diagm(ones(size(u,1)))),inds=1:size(u,2))
# 	nRC = size(u,1)
# 	nInd = length(inds)
# 	df,Ψ0 = params(prior)
# 	Ψ = Matrix{Float64}(Ψ0)
# 	eps_i = zeros(nRC)
# 	for ii = inds
# 		@views eps_i .= u[:,ii] #if mean is not zero, would need to subtract mean here!
# 		mul!(Ψ,eps_i,eps_i',1.0,1.0) #Ψ += eps_i*eps_i'
# 	end
# 	posterior = InverseWishart(I+df,Ψ)
# end

########################################################
## update utilities and availability
########################################################
mylog(x) = x > 0 ? log(x) : -Inf

function updateIndividual!(theta,lv,cc,temps,ii)
	myid = Threads.threadid()
	rng = myRNG[myid]
	J = size(lv.U,1)
	Ui = view(lv.U,:,ii)
	Ai = view(lv.Astar,:,ii)
	feasi = view(cc.feasibleUnranked,:,ii)
	u0 = view(lv.U0,:,ii)
	enrolli = cc.enroll[ii]
	typei = cc.typ[ii]
	#human-capital index H
	indEnroll = temps.indEnroll[ii]
	if enrolli > 0 && !cc.specialAdmit[ii]
		muH = dot(view(temps.Xoutcome[typei],indEnroll,:),view(theta.betaOutcome,:,typei))
		if cc.Grad6[ii]
			lv.H[ii] = rand(myRNG[myid],truncated(Normal(muH,1.0),0,Inf))
		else
			lv.H[ii] = rand(myRNG[myid],truncated(Normal(muH,1.0),-Inf,0))
		end
	else
		lv.H[ii] = NaN
	end
	#= outside options
	At application time, oo is u0[1].
	At enrollment time, oo is max(u0[1],u0[2]).
	Because enrollment-time constraint is on max, we have:
		(a) neither outside option can exceed aftermarket upper bound oo_ub_a
		(b) at least one outside option must beat aftermarket lower bound oo_lb_a.
	Point (b) means that u0[k] has to beat oo_lb_a only if u0[-k] does not.
	app constraint applies only to first outside option. =#
	oo_lb_app = getApplicationBound(Ui, feasi, u0, cc.OO_lb[ii],false,true)
	oo_ub_app = getApplicationBound(Ui, feasi, u0, cc.OO_ub[ii],true,true)
	(oo_lb_a, oo_ub_a) = getAftermarketBounds(Ui,u0,Ai,0,enrolli) #bounds max of u0[1], u0[2]
	oo_lb = max(oo_lb_app,oo_lb_a)
	oo_ub = min(oo_ub_app,oo_ub_a) #UB: must satisfy both upper bounds.
	lb0 = (u0[2] > oo_lb_a) ? oo_lb_app : oo_lb #lower bound on first oo depends on 2nd oo
	lb0 < u0[1] < oo_ub || error("bounds cross in market $(cc.year), ii=$ii, u0[1]. (lb=$lb0, u0[1]=$(u0[1]),ub=$oo_ub)")
	u0[1] = rand(myRNG[myid],truncated( Normal(temps.mu0[1,ii],sqrt(theta.sigsqOO0[typei])), lb0, oo_ub))
	lb1 = (u0[1] > oo_lb_a) ? -Inf : oo_lb_a
	u0[2] = rand(myRNG[myid],truncated( Normal(temps.mu0[2,ii],sqrt(theta.sigsqOO1[typei])), lb1, oo_ub_a))
	#now do utility and availability of inside options
	for jj=1:J
		lb_app = getApplicationBound(Ui, feasi, u0, cc.J_lb[jj,ii],false)
		ub_app = getApplicationBound(Ui, feasi, u0, cc.J_ub[jj,ii],true)
		(lb_a, ub_a) = getAftermarketBounds(Ui,u0,Ai,jj,enrolli) #bounds on utility if program were available
		if cc.waitlistMightCall[jj,ii] #formerly, cc.A_ub[jj,ii] > 0
			alpha = view(theta.alpha,:,typei) #cc.platform[jj] ? theta.alphaOn : theta.alphaOff
			xf = view(cc.Xfriction_full[ii],jj,:)
			Ai[jj] = drawAvailability!(Ui[jj],lb_a,ub_a,cc.A_lb[jj,ii],cc.A_ub[jj,ii],xf,alpha)
		end
		if Ai[jj] > 0
			lb = max(lb_app,lb_a)
			ub = min(ub_app,ub_a)
		else
			(lb,ub) = (lb_app,ub_app)
		end
		lb < Ui[jj] < ub || error("bounds cross in market $(cc.year), ii=$ii, jj=$jj: Ui[$jj]=$(Ui[jj]), (lb,ub)_app=($lb_app,$ub_app), (lb,ub)_after=($lb_a,$ub_a)")
		if enrolli == jj && !cc.specialAdmit[ii]
			Ui[jj] = draw_u_enrolled!(theta,lv,temps,jj,ii,typei,lb,ub)
			temps.Xoutcome[typei][indEnroll,1] = Ui[jj]
		else
			Ui[jj] = rand(rng,truncated(Normal(temps.mu[jj,ii],1.0),lb,ub))
		end
	end
end

# function draw_u0_enrolled!(theta,lv,temps,ii,typei,lb,ub)
#     a = theta.betaOutcome[1,typei]
# 	Xoutcome = temps.Xoutcome[typei]
# 	ind = temps.indEnroll[ii]
#     Xoutcome[ind,1] = 0
#     muH = dot(view(Xoutcome,ind,:),view(theta.betaOutcome,:,typei))
#     htilde = lv.H[ii] - muH
# 	sig2tilde = 1/(a^2 + 1/theta.sigsqOO0[typei])
#     mutilde = (temps.mu0[1,ii]/theta.sigsqOO0[typei] + a*htilde)*sig2tilde
#     uij = rand(myRNG[Threads.threadid()],truncated(Normal(mutilde,sqrt(sig2tilde)),lb,ub))
# end

function draw_u_enrolled!(theta,lv,temps,jj,ii,typei,lb,ub)
    a = theta.betaOutcome[1,typei]
	Xoutcome = temps.Xoutcome[typei]
	ind = temps.indEnroll[ii]
    Xoutcome[ind,1] = 0
    muH = dot(view(Xoutcome,ind,:),view(theta.betaOutcome,:,typei))
    htilde = lv.H[ii] - muH
	sig2tilde = 1/(a^2 + 1)
    mutilde = (temps.mu[jj,ii] + a*htilde)*sig2tilde
    uij = rand(myRNG[Threads.threadid()],truncated(Normal(mutilde,sqrt(sig2tilde)),lb,ub))
end

function drawAvailability!(u,lb,ub,alb,aub,xf,alpha)
	# alb==aub && return alb
	# lb < u < ub || return false
	lowerbound = alb==0 ? -Inf : 0.0
	upperbound = (aub==0 || u > ub) ? 0.0 : Inf
	#rand(myRNG[Threads.threadid()]) > alpha #alpha is probability of UNavailability
	rand(myRNG[Threads.threadid()],truncated(Normal(dot(xf,alpha),1.0),lowerbound,upperbound))
end

function getAftermarketBounds(Ui,u0,Ai,jj,enrolli)
	u0_post = max(u0[1],u0[2])
	if enrolli > 0 && jj != enrolli
		ub = Ui[enrolli]
		lb = -Inf
	elseif enrolli == 0 && jj != enrolli
		ub = u0_post
		lb = -Inf
	elseif enrolli == jj
		ub = Inf
		lb = -Inf
		for kk=1:length(Ui)
			(Ai[kk]>0) && (kk != enrolli) && (lb = max(lb,Ui[kk]))
			enrolli > 0 && (lb = max(lb,u0_post))
		end
	end
	(lb,ub)
end

function getApplicationBound(Ui,feasi,u0,k,upper=true,isOO=false)
	if k>=1
		return Ui[k]
	elseif k==0
		return u0[1]
	elseif k==-1 || (k==-2 && !any(feasi))
		return ifelse(upper,Inf,-Inf)
	elseif k==-2
		bound = isOO ? -Inf : u0[1]
		for kk=1:length(Ui)
			feasi[kk] && (bound = max(bound,Ui[kk]))
		end
		return bound
	end
end

function updatemu!(theta,lv,cc,temps,ii,programRE=true,fixed=true,ij=true,rc=true)
	mui = view(temps.mu,:,ii)
	typei = cc.typ[ii]
	mui .= 0
	programRE && (mui .+= view(temps.ure,:,typei)) #view(theta.programRE,:,typei))
	fixed && (mui .+= view(temps.ufe,:,typei))
	ij && (mui .+= view(temps.umatch,:,ii))
	rc && (mui .+= view(temps.urc,:,ii))
	nothing
end

function justUpdateIndividuals!(theta,dataset,priors,T)
    for year in keys(dataset)
        (cc,lv,temps) = dataset[year]
        precomputeMeanUtilities!(theta,lv,cc,temps)
    end
    for t=1:T
        for year in keys(dataset)
            (cc,lv,temps) = dataset[year]
            (J,I) = size(lv.U)
            println("year $year, utility update $t")
            Threads.@threads for ii = 1:I
                updatemu!(theta,lv,cc,temps,ii)
                updateIndividual!(theta,lv,cc,temps,ii) #(1) utility and availability
                updateRC!(theta,lv,cc,temps,ii) #(2) students' random coefficients on program char.
            end
        end
    end
end
